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We present a study of spinning black hole binaries focusing on the spin dynamics of the individual 
black holes as well as on the gravitational recoil acquired by the black hole produced by the merger. 
. . . We consider two series of initial spin orientations away from the binary orbital plane. In one of 

' the series, the spins are anti-aligned; for the second series, one of the spins points away from the 

\ binary along the line separating the black holes. We find a remarkable agreement between the spin 

. dynamics predicted at 2nd post-Newtonian order and those from numerical relativity. For each 

' configuration, we compute the kick of the final black hole. We use the kick estimates from the series 

^ with anti-aligned spins to fit the parameters in the Kidder kick formula, and verify that the recoil 

^ , along the direction of the orbital angular momentum is oc smO and on the orbital plane oc cos0, 

' with 9 the angle between the spin directions and the orbital angular momentum. We also find that 

the black hole spins can be well estimated by evaluating the isolated horizon spin on spheres of 
constant coordinate radius. 
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bJO, I. INTRODUCTION 

CN ' Immediately after the discovery of the Moving Puncture Recipe (MPR) [3, [2| , a recipe providing the ingredients to 
' successfully evolve binary black holes (BBHs), the numerical relativity efforts focused on studying the gravitational 
' . ^ ■ recoil or kick acquired by the black hole (BH) produced in the merger 0, 0, IB| ■ The main driving force behind these 
^ studies has been the astrophysical implications of these kicks on the supermassive black holes (SMBHs) at the centers 
, of galaxies [1,01 • Specifically, a detailed understanding of these kicks is vital to explain the demo grap hics, growth 
and merger rates of SMBHs H, Q, as well as their absence in dwarf galaxies and stellar clusters p^.llll|. 

When viewed in terms of modes of the gravitational radiation emitted by the binary, kicks arise from the overlap 
of those modes [l^ [l3| . A non- vanishing overlap will be produced if the BHs in the binary have un-equal masses 
', and/or are spinning with non-trivial relative orientations. For kicks from non-spinning BBHs, the most comprehensive 
numerical relativity study 0] showed that one can parameterize the magnitude of the kick velocity as 



O 



H (1 + qY 



l + B- 



(1) 



with A = 1.2 X 10^ km s"""^ , B = —0.93 and q — -M1/M2. This parameterization was motivated by the scalings 
originally introduced by Fitchett [3, From Eq. ([1]), the maximum kick has a magnitude of 175 km s~^ and 
occurs Sit q — 0.36 or symmetrized reduced mass 77 = Mi M2/M^ = 0.195, with M — Mi -f M2 the total mass of the 
binary. Other mass parameters that will be used are 5M = Mi — M2 and = Mi M2/M. When compared to the 
escape velocities of galactic structures, the kicks from non-spinning and un-equal mass binaries are modest. They are 
not high enough to eject the BH from its host galaxy (Tlj . 

The next frontier was to investi gate kicks in which the emission of linear momentum was due to the spin of the 
BHs. The first study of this kind [lj| produced kick velocities of V — 475 km s^^ a/M for BHs with opposite and 
equal magnitude spins parallel to the orbital angular momentum. Similar studies followed soon after ^Q, [13] that 
produced complementary results. The kick velocities of ~ 500 km obtained from these configurations could in 
principle explain the absence of massive BHs in dwarf ellipticals [ll| . Motivated by post-Newtonian (FN) results [3 1 
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it was immediately realized that the orientation of the BHs' spins has a profound effect on the kick that the final BH 
receives. Gonzalez et al. [l^ carried out the first simulations in which the spins of the BHs are initially anti-aligned 
in the orbital plane and found that kick velocities of at least 2500 km are possible. Similar studies [l^l suggest 
that the kick could be scaled to reach a maximum of ~ 4000 km s~^ . 

As more studies of gravitational recoil continued to emerge, generalizations of the phcnomcnological kick formula 
Eq. ([1]) to include spins have been introduced [13, HH, , all motivated by the structure of the formula for the rate 
of linear momentum radiated, a formula first derived by Kidder p^ . The terms involving spin-orbit effects in this 
formula read 

^ = 'l5~JlTW ^ ^) ~ (n X S) - (n X v) [3r(nl]) + 2(^;S)]} , (2) 

where (ab) denotes the vector dot product, i.e. (ab) = a • b. We are following as closely as possible the notation in 
Ref. [23| and introduce the spin variables 

S = S1+S2 

S2 Si 



where the vector x denotes the relative position vector of M2 with respect to Mi, with r = |x|, v = dx/dt, n = x/r 
and Ln = /i x x v, the Newtonian angular momentum. We also introduce a flat-space orthonormal rotating triad 
{n, k, 1} such that k = 1 x n with 1 = Ln/|Ln| and hence 1 is perpendicular to the orbital plane. 
With these definitions, Eq. ([2]) has the following structure: 

^ = [...](kxS) + [...](nxS) + {[...](fcE) + [...](nS)}l, (3) 

or equivalently 

^ = {[...]k+[...]n}(/S) + {[...](fcE) + [...](nS)}l, (4) 

where we only show the explicit dependence on S relative to the orthonormal tetrad. Given the form of Eq. ([4]), we 
propose the following parameterization of the contribution of the spins to the gravitational recoil: 



V = 



M2 + 



-[[Hkk + Hnn]ila) + [Kkika) + Knina)]l'j , (5) 



where a = S/|S|. We will refer to Eq. ([5]) as the Kidder kick formula.^ The parameters Hk, Hn, Kk and Kn in 
Eq. (O are to be determined from numerical simulations. A fundamental aspect of the validity of this formula is the 
dependence of the kick velocity on the cosine angles {Icr), (ka) and (na). Spin precession will force these angles to 
evolve in time. Thus, one is faced with the task of measuring the entrance angles. These are the angle values when 
the binary reaches the "last" orbit or plunge, namely the time that signals the beginning of the phase when the bulk 
of the kick gets accumulated. An identification of the entrance angles would allow one to determine the Hk, Hn, Kk 
and Kn parameters in Eq. ([5|) from numerical simulations. 

The work in this paper is aimed at exploring the parameter space of spinning BBHs with focus on the dynamics of 
the individual spins and the kick that the final BH receives. We consider two series of equal mass BHs (i.e. SM = 0). 
In one series, called the B-series, the BHs initially have equal spin magnitudes and anti-aligned directions. That is, 
S = and S = 4S2 = —4 Si. The elements of this series arc obtained by changing the orientation of S relative to 
the unit vector 1. In the second series, called the S-series, we also keep the spin magnitudes constant. What changes 
in this series is the relative alignment of the spins. For each run in both series, we monitor the precession dynamics of 
the individual spins and compare them with PN predictions. We find a remarkable agreement with 2PN results: the 
2PN dynamics closely match those from numerical relativity up to the point when a common apparent horizon (AH) 
is formed. For all models, we compute the gravitational recoil on the final BH. We use the kick estimates from the 
B-series to find parameters in the Kidder kick formula and also verify the angular dependence in V that this formula 



^ There are several versions of parameterized kick formulas. Since all are motivated by Kidder's seminal work [181] , we will generically call 
them Kidder kick formulae. 
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implies. As numerical relativity efforts explore different regions of the parameter space, the values of the parameters 
in Eq. ([5]) will be improved or validated. A phenomenological formula of this kind is of great value for astrophysical 
studies such as those explaining the population of SMBHs. 

The paper is organized as follows: In Sec. |TT1 we use a multipole analysis to demonstrate the dependence of the 
kicks on the spin orientations as given by the Kidder kick formula. In Sec. IIIIl we summarize our computational 
infrastructure. A detailed description of the two series of initial data configurations is given in Sec. IIVI The analysis 
of the BH spin dynamics is presented in Sec. |Vl Kick results, including the fit to the Kidder kick formula, are given 
in Sec. |Vll We end with conclusions in Sec. IVIII 



II. KICKS AND ENTRANCE ANGLES 



To gain further understanding of the Kidder kick formula, we present an analysis based on the multipole formulas 
of Refs. [13 1 in which the rate of radiated linear momentum is estimated, to lowest order, as an interference of 
the mass and spin quadrupoles. Excluding non-spin terms, this formula reads 

^ _ ]^.jk r(3) „(3) , ± „(4) „(3) J_ «jfe .(4) rrii) /rN 
- ^jl ^kl ^ Q^^ijk^jk 12Q ^Jlm^klrn- 1°^ 

Here lij and lijk are respectively the mass quadrupole and octupole. Similarly, Hij and Hijk are the spin quadrupole 
and octupole, respectively. In Eq. a super- index ("^ denotes an nth-time derivative. 

In previous work , we used the first term (interference between the mass and the spin quadrupoles) to estimate 
the kick from quasi-circular inspiral to merger by integrating Eq. ([G]). This term is periodic, with period equal to 
the orbital period, so the kick is dominated by the "last" half orbit in the inspiral. The estimate is computed by 
integrating over a close-in half orbit (as in Section |T1 the result depends on the magnitude and direction of the spins 
with respect to the orbital angular momentum L = LI), and absorbing the resulting error as a normalization constant, 
where the constant is fixed by comparing estimate to numerics for one configuration. We take the same approach 
here. 

Note that the second term in Eq. ^ will be quadratic in the spin, but the spin multipoles have one extra factor 
(81^2/^/1, 2)/^ (where d is the "last orbit separation", and of order several M) that suppresses the radiation from this 
term by the same factor compared to the first term. While this term's contribution may become important in the 
future, for the moderate spin values we (and others) are currently considering, we do not expect significant nonlinear 
dependence. The third term vanishes (the mass octupole vanishes) for equal mass circular orbits as appropriate to 
our computational quasi-circular inspiral, so the equation in our current context is just the first term. 

For the purpose of investigating the entrance angles, we consider a binary system consisting of equal mass BHs in 
circular orbit initially confined to the xy plane. The orbit is initially oriented so that the BHs are located on the 
X-axis, the BHi on the positive cc-axis and BH2 on the negative a;-axis. We discuss first the case in which only the 
BHi is spinning. We parameterize the orientation of the spin using the usual (fixed frame) polar and axial angles 9 
and tp. Thus we have Six = Si smOcos(p, Siy — Si sin sin and Siz = Si cosO. 

The calculation of the mass quadrupole is straightforward, see e.g. [l^; the spin quadrupole can be most easily 
calculated by imagining a spin dipole (charges ±Mi/2, separation Si/Mi) and conceptually taking the hmit at the 
end. The spin enters only linearly in the spin quadrupole Hki- The structure is different for spin components in 
different directions, and we can compute them independently for the different components. The nonzero components 
are: 

For Six-. 

^^^Hf^ = ^-dSixUo'siniu^t) 
^""^H^f = -^d Six u;"" cos {cut) 
^"^^yy = -IdSixUJ^'smiiot) 

^^^Hil'> = -^dSixUJ^'smiujt); (7) 
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For Si, 



^"^H^Ix = Sly cos {ujt) 
= -]-dSiyW^sin{ujt) 



For Si, 



3 

(^^iJi^) = id5i^c^3^osM; (8) 
^"^^^^^^ = ]^dSi,Lo^sin{ujt) 

^'^H^f - -id5i,c^3^osM). (9) 

The spin quadrupole for arbitrary spin direction is the sum of the Six, Siy, Siz terms. In deriving these expressions, 
we assume that spins, which are parallel transported in the evolution, remain constant in Cartesian coordinates. 
This approximation is adequate for the level of accuracy of these estimates. The radiated linear momentum equation 
Eq. © is then explicitly: 

J pa; Q 

^ = ^M^d^uj^Si.siniiot) 

^ = — ^M^d^^e^i.cosM) 
at 45 

—— = -—M^d^Lu^\SixCOs(ujt)~SiySm(Lut)]. (10) 
45 

The in-plane component of the force rotates with the orbit; the out of plane component oscillates at the frequency of 
the orbit. 

If there is a spin on the second hole, the forms are the same as Eqs. (fTU)) . but the angle ujt is replaced by uut + tt. 
This replacement has the effect of introducing a global minus sign into the spin quadrupole for the spin on the hole 
initially located on the negative a;-axis. This means that the kick estimate is doubled if the second spin is equal and 
opposite, but we estimate zero kick if the spins are equal and parallel. For generic second spin, as in our S-series, one 
simply subtracts the components in Eqs. (|10p for this second spin 5*2 from those for the first. 

We concentrate our attention on the B-series. As we shall see later, this is the series for which we are going to 
be able to verify, from our simulations, the dependence of the Kidder kick formula on the entrance angles. In the 
B-series, the spins are fixed magnitude. Hence, the Eqs. (|10p . including the contributions from both spins, read: 

^ ^ — M^d^u^ Si cose sin (ut) 
dt ^ ' 

= M'^d^ uo^ Si cose cos (tut) 

dt A5 ^ ' 

dP^ 32 

— — = -— M^d^uj^ Si sine cos iujt + ip). (11) 
dt 45 

Eqs. pT|) predict a z-kick V^- cx sin0 and kicks oc cos0 in the orbital plane. Notice also the dependence of the z-kick 
on the entry angle [ut-'rip), demonstrating the fact that the net z-kick can vanish for carefully chosen entry angle. For 
the circular orbits treated here, dP^ /dt in Eqs. p4|) identifies the quantities Kk and Kn in the Kidder kick formula, 
Eq. (O, as equal. We will compare the predictions of Eqs. pT|) on the scaling of the kicks with the angle e in Sec.lVIl 



III. COMPUTATIONAL METHODOLOGY 



We follow the MPR to evolve the BBH configurations. Briefly, the MPR builds upon the BSSN system of evolution 
equations [13, HI, HI] , models BHs with "punctures" [13] and uses dynamic gauge conditions [3, 01 designed to allow 
these punctures to move. The explicit form of the evolution equations for the lapse and shift gauge quantities are the 
"covariant" form of the "l+log" slicing {dt — P^di)a = —2aK and a modified gamma-freezing condition [29l [30j 
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for the shift: 9t/3* = B^,B^ = dtT^ —S.dtP'' — P^djV^, where K is the trace of the extrinsic curvature, the trace of the 
conformal connection and ^ = 2 a free, dissipative parameter. The importance of these gauge conditions is twofold: 
First, they avoid the need of excising the BH singularity from the computational domain since they effectively halt 
the evolution (i.e. lapse function a vanishes) near the BH singularity 'si'l. Second, they allow for movement of the BH 
or puncture through the computational domain while freezing the evolution inside of the BH horizon. See Ref. (3^ 
for a detailed description and analysis of the MPR. 

Our source code was produced by the Kranc code generation package [33l| and uses the Cactus infrastructure [s^l 
for parallelization and Carpet |35j for mesh refinement. The code uses fourth order accurate finite differencing 
(centered for all non-advection and a lop-sided stencil for the advection terms) and a fourth order Runge-Kutta 
temporal updating scheme with Courant factor of 0.5. The initial data code was developed by Ansorg et al. [3y|. The 
initial free parameters (e.g . sp ecifying angular momentum, spins, masses, separations) are chosen according to the 
effective potential method [37l . Issj or usingPN parameters [3, [33] ■ These methods both yield BBH initial data sets 
representing BBHs in quasi-circular orbit |32l |. 

The computational grids consist of a nested set of 10 refinement levels, with the finest mesh having resolution 
h = M/35.2. This resolution translates into a resolution of about h = to/14, with respect to the bare mass, to, of the 
punctures according to Tables [J and [ill The minimal resolution found to be adequate for spinning cases according 
to Campanelli et al. [3^ is h < M/30. The grid sizes in our h — M/35.2 simulations are: the 4 finest refinement 
levels have 44^^ grid-points plus 6 coarser refinement levels of 88"^. All grids are initially cubical. During the evolution, 
the shape and number of grid-points per refinement level vary due to adaptivity. The coarsest mesh is kept fixed and 
extends to 640 M from the origin in each direction. Because the simulations in this work are very similar (regarding 
mesh setups, grid sizes and refinement scales) to those in our previous work fisj, the convergence and errors estimates 
in the present study are comparable. 

In order to study the spin dynamics of the BHs, we need infrastructure to compute the individual spins of the BHs. 
The isolated horizon formalism [i^ [4l|, \^ provides a definition associated with a Killing vector of the spacetime of 
the spin of a single BH: 

Sv = Z- { v'n'K,jdS (12) 

iSTT J AH 

where if^ is a Killing vector on the AH surface, Kij is the extrinsic curvature of the 3D-slice and is the outward 
pointing unit normal vector to the AH. The direction of the spin is given by the Killing vector ip^ . To facilitate finding 
the spin direction, Campanelli et al. 43] introduced the usage of the fiat space coordinate rotational Killing vectors 

(0,-5,y) 
(i,G, -x) 

where the coordinates {x,y,z) are relative to the position of the BH. The spin is then given by S = (Sx, Sy, Sz), 
where each component is obtained, in the fixed {x, y, z} coordinate system, by evaluating Eq. ()12l) with each of the 
coordinate rotational Killing vectors. There is an excellent agreement between the approximate spin this method 
yields and the one using the Killing vector (/p* (when one exists) [isf. There are efficient AH finders |4J] available; 
however, they impose a non-negligible overhead in the simulations. To gain efficiency, we relax the condition that the 
integral in Eq. has to be evaluated at the AH and choose a coordinate sphere around the puncture. The radius 
of the sphere is chosen sufficiently small, that the sphere is contained within the BH's horizon. 

Fig. [B shows a comparison of the Sx component between the values using the AH surface and three different 
coordinate spheres of radius r for the S-90 model (see Table [Hi BBH evolution. There is good agreement into the 
merger regime. The vertical line in Fig. [T] and subsequent figures shows the first time a common AH is found. 
After that time, no individual apparent horizons exist and the spheres centered on the punctures track different and 
meaningless values of Sx ■ 



IV. INITIAL DATA AND RADIATED QUANTITIES 

We consider two series of equal mass BHs (i.e. 5M = 0). In both series, initially the BHs have the same spin 
magnitude Si/Mf = S2/M2 = 0.6. The initial orientation of the BH's spin is in the a;z-plane. We use the polar angle 
61: the angle between the z-axis and the direction of Si in the xz-plane. The BH2 is treated similarly. We take the 
convention that positive (negative) 9 angles are measured (counter-) clockwise in the xz-plane. In the series referred 
to as B-series, the BH's are anti-aligned, i.e. 9 = 02 ^ &!- 180°, so Si = -S2. That is, S = and S = 4 S2 = -4 Si. 
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time (M) 



FIG. 1; Comparison of computed on the horizon and from spheres with radius r for a BBH evolution (S-90 model, see 
Table inil. The vertical line (h ere and in subsequent Figures) shows the first time a common AH is found. 

model x[M] py[M] y[kms"^] J,,d[%I^g] grad[%M] r.„,.[M] 



B-20 


2.986 


0.138 


427 


24 


3.3 


109.1 


B-30 


2.990 


0.138 


544 


24 


3.3 


109.1 


B-50 


3.000 


0.137 


761 


25 


3.4 


108.6 


B-70 


3.009 


0.137 


908 


25 


3.4 


108.6 


B-80 


3.012 


0.137 


945 


25 


3.4 


108.4 


B-90 


3.013 


0.137 


963 


25 


3.4 


108.4 



TABLE I: B-series: Initial data parameters for the B-series. The models in this series are labeled as B-0, where the angle 
d = 02 ~ 9i — 180° {9 — 0° corresponds to spins parallel and anti-parallel to the orbital angular momentum). The punctures 
have bare masses mi, 2 ~ 0.395, are located on the x-axis at and have initial momentum ^^py in the ^/-direction. Results 
listed are the magnitude of the recoil velocity V, the radiated angular momentum J^^^ in % of the initial orbital angular 
momentum LJ, the energy radiated -Erad, and the time Tmax which is an estimate of the merger time derived from the time it 
takes in each simulation to reach the maximum amplitude in 'if 4. 

The elements in this series are obtained by changing 9. In the S-series, we initiaUy orient Si to 9i = 270° ~ —90° 
and vary 9 = 92- 

We chose orbital parameters (i.e. bare masses, separation and momentum) in the B-series by minimizing the 
effective binding energy [H, H^, while for the S-series we used PN parameters [l^ H^. Initially, BHi is located 
at position (— x/M, 0,0) and has linear momentum (0, — Pj,/Af , 0). Similarly, BII2 is at position (x/M, 0,0) with 
linear momentum (0,Pj,/M, 0) . It turns out that the bare puncture masses for both series are roughly constant, 
TOi = m2 « 0.395 M to the 3rd digit of precision. The slight changes are needed to keep the irreducible masses Mi = 
M2 = 0.5 AI. As mentioned above, the spins in both BHs are initially in the xz-plane; that is, (5f 0, Sf 

where S'f 2 = 5'i,2 sin6'i,2 and S'f 2 = 'S'1,2 cos0i,2 with 51,2 = 0.15 M^. 

Table Jllists the relevant initial data parameters for the B-series, while Table HIl gives the parameters for the S-series. 
In addition to the initial data parameters, the tables also report the radiated angular momentum J^^^ in % of the 
initial orbital angular momentum, L^, as well as a time estimate of the common AH formation. We use the maximum 
in ^'4 shifted by the extraction radius and an additional 10 Af as an indicator for the merger time Tmax- We have 
found that this measure is accurate to a few M. For the B-series, the spin of the final BH is J/AP — 0.62 for all 
models. Constant in both series is the total ADM mass, Eabm ~ 0.985 M. While the runs B-90 and S-90 have the 
same spin configurations, i.e. spins pointing along the x-axis only, the radiated energy and angular momentum are 
different because they differ in initial separation and angular momentum. The radiated quantities were extracted 
at r = 40 Af. For a number of models, we have carried out simulations at lower resolution (A//32) and measured 
at detector radii r/M = {30,40,50,60,80}. Based on the variations observed in the measured quantities (energy, 
angular momentum and kicks), we estimate the reported numbers to be accurate to about 15%. 
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model mi,2[M] Py[M] \/[km ] J,ad[%Lg] grad[%M] Jg„,i[M^] Tn.ax[M] 



S-0 


0.396 





132 


854 


34 


4.6 





68 


192.3 


S-15 


0.396 





132 


1401 


33 


4.4 





68 


189.5 


S-30 


0.396 





132 


2000 


33 


4.4 





67 


184.1 


S-45 


0.396 





133 


2030 


32 


4.3 





66 


177.3 


S-60 


0.395 





134 


1218 


30 


4.0 





65 


168.6 


S-75 


0.395 





135 


230 


28 


3.7 





64 


159.1 


S-90 


0.395 





137 


1462 


26 


3.4 





62 


148.6 


S-105 


0.395 





138 


1979 


25 


3.3 





60 


138.6 


S-120 


0.395 





139 


1787 


24 


3.2 





58 


130.5 


S-135 


0.395 





140 


1234 


23 


3.0 





56 


124.1 


S-150 


0.395 





141 


689 


21 


2.9 





55 


119.5 


S-165 


0.395 





141 


335 


21 


2.8 





55 


117.7 


S-180 


0.395 





141 


188 


20 


2.8 





55 


117.7 


S-195 


0.395 





141 


157 


20 


2.8 





55 


120.5 


S-210 


0.395 





141 


173 


22 


3.0 





56 


125.5 


S-225 


0.395 





140 


223 


22 


3.2 





57 


132.7 


S-240 


0.395 





139 


268 


23 


3.4 





59 


141.4 


S-285 


0.395 





135 


253 


26 


3.9 





65 


174.1 


S-300 


0.396 





134 


406 


29 


4.2 





66 


181.8 


S-315 


0.396 





133 


399 


31 


4.5 





67 


187.7 


S-330 


0.396 





132 


354 


32 


4.6 





68 


191.8 


S-345 


0.396 





132 


459 


33 


4.6 





68 


193.2 



TABLE 11: The S-series. For all cases, initially the BHi is located along the a::-axis at a:: = —3.1 M, has momentum pointing 
along the y-direction with value — Py, and has spin Si = (—0.15 /A/'^,0,0), thus 9i — —90° and ipi = —180°. BH2 is located 
also along the a;-axis but at x = 3.1 M with momentum py. In these runs, labeled S-9, the angle 6 gives the angle in the 
a;z-plane that the spin of BH2 makes with respect to the z-axis. Results listed are the magnitude of the recoil velocity V , the 
radiated angular momentum J^^^ in % of the initial orbital angular momentum LJ, the energy radiated i?rad, the spin of the 
final BH Jlnai along the z-axis, and the time Tmax which is an estimate of the merger time derived from the time it takes in 
each simulation to reach the maximum amplitude in 'I'4. 



SPIN DYNAMICS 



In the present work, we are interested investigating the degree to w^hich the spin dynamics described by PN equations 
agrees with that from numerical relativity. Following Ref. [46l] . the precession equation of BHi in the binary with 
mass Afi, spin Si, position Xi and velocity vi is given by 



at 



(13) 



which imphes that BHi precesses around the vector ili with rate The precession angular frequency vector fii 

is given to 2PN by 



Oi = 



M2 
^2 



-ni2 X vi — 2ni2 x V2 



, „, ^2 1 2 / ^ 2 7 Ml IM2 
ni2 X Vi - - ("12^2 ) + o ^1 ^ (^1^2) + ^2 + 7: 

9M2 
2~V 



ni2 X V2 3(^12^2)^ + 2{viV2) -2vl + + 



+ vi X V2(^3(ni2Wi) - ^(^12^2)) 



(14) 



with X = Xi — X2, r = |x| and ni2 = x/r. The expressions for the companion BH2 are obtained by switching 1 2 in 
Eqs. (|13m4p . In Eq. (|14p . the first term in square brackets represents the IPN contribution. For comparison, we also 
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110 120 130 140 150 

time (M) 



FIG. 2: Comparison of dSi/dt computed from the numerical evolution directly and by using PN formulas for the S-90 run. 
Kidder describes the dynamics using precession angular frequency given by Eq. (I15|) . Blanchet IPN denotes the dynamics with 
f2i given by the first term in Eq. (|14p: Blanchet 2PN denotes the case in which the entire expression in Eq. (|14|l is used. The 
vertical line around t = 149M indicates the formation of a common apparent horizon. 



show the precession angular frequency from Kidder |18l |. in which the terms oc L^r (corresponding to the first line in 
Eq. are accurate to IPN but the expression also contains spin-spin terms: 



^1 = ^ 



-'N 



2 Ml' 



3(71i2S'2)ni2 



(15) 



where Lat = jix x Vi2 denotes the Newtonian angular momentum. Here again one obtains the expression for BH2 by 
switching 1^2. 

Fig. [5] shows the time evolution of dSi/dt computed in four different ways for the S-90 run (the spins are equal 
magnitude and anti-aligned in the xy-plane). The time evolution of 6,^2 /dt is equal and opposite in this case. Solid 
lines, labeled numeric, represent the numerical relativity solutions. The values of dSi/dt are obtained by constructing 
each BH spin as described in Sec. IIIIl followed by finite differences to approximate the time derivative. A long dashed 
line, labeled Kidder, denotes dSi/dt computed using the precession angular frequency Eq. p^ . The dotted line, 
labeled Blanchet IPN, represents the result from using only the IPN contribution in the precession angular frequency 
Eq. (HU; that is, it corresponds to Kidder's precession without the inclusion of spin-spin interactions. Finally, the 
dashed-dotted line, labeled Blanchet 2PN, depicts the evolution oidSi/dt using the entire expression in Eq. ([14]). In 
the construction of the PN precession angular frequencies, we use the positions and velocities of the punctures from 
the numerical simulations. The vertical lines in Fig. [2] denote the time at which a common AH is formed. 

It is remarkable how accurately the 2PN approximations of dSi/dt track the numerical result deep into the merger 
regime, close to the formation of a common AH. Comparisons beyond the time when a common AH forms are not 
very meaningful since the individual trapped surfaces loose their horizon interpretation and our spin measure breaks 
down (see Sec. IIIIl) . Also interesting is that the spin-spin terms in Kidder's expression make only a small contribution 
to dS^ /dt and dS^ /dt, as can be seen from the similarities of the Kidder and Blanchet IPN lines. On the other 
hand, the spin-spin are responsible for the differences between the Kidder and Blanchet IPN values of dS^ /dt near 
the mergers, as one can observe in the bottom panel of Fig. [21 This discrepancy can be traced to the 2;-component in 
the third term in Eq. (jlSp . The first term when using the frequency Eq. (jlSp in Eq. (|13p contains the z-component 
of L^r X Si, which is numerically very close to zero for the S-90 model. The second term contains the z-component 
of S2 X Si, which is also close to zero. In the third term, we have {ni2S2) and the z-component of ni2 x Si. Both of 
these terms grow rapidly near the merger. In particular, the z-component of ni2 x Si develops significant noise which 
terminates the line early. Finally, it is very clear that including terms up to 2PN makes an important difference in 
improving the matching to the numerical solution. 

We have carried out comparisons similar to that in Fig. [2] for all the runs in the B- and S-series. The results in 
every case are similar; namely, that the dynamics of dSi^2/dt are very well approximated by 2PN, and this description 
only starts breaking down close to the merger. 

The most significant variation observed in the dSi^2/dt dynamics from one S-series run to another was the time 
at which the merger takes place or, equivalently, the time at which the gravitational radiation emitted reaches its 
maximum (see Table [ll]). The differences in merger time are due to spin hangup 43] of the merger. Figs. [3|[7] show 
the comparison for a few selected models from the S-series. The plots show dS/dt for both BHs. The left panel shows 
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FIG. 3: Comparison of S-15 run numerical to Blanchet 2PN. Left panel shows the results of the comparison for BHi and the 
right panel for BH2. The top plots on each panel show with a solid line dS^/dt from our numerical simulations and with a dashed 
line the values from Blanchet 2PN. The labels denote each component. The bottom plots on each panel show the difference 
between the numerical solution and the Blanchet 2PN, with solid, dashed and dotted lines for the x, y and z components, 
respectively. 

BHi, which has the initial spin direction along the negative a;-axis. In this case, there is good agreement between the 
2PN and computational results for all models. In the right panel, we show BH2 with a different initial spin direction 
specified according to the angle 9. For spins of BII2 more parallel to the orbital angular momentum, the precession 
becomes smaller (note range of the y-axis to the right of the figure). It seems very likely that with the small precession 
shown in the S-15 model, the visible disagreement to the PN result is just a numerical artifact that could be cured 
by higher resolution. 

To further understand the spin dynamics, we focus our attention to the evolution of the direction of the spins Si, 2 
and the vector S. Fig. [5] shows the evolution of the spin directional angles 61^2 and 931,2 for the B-series. The angles 
and ip are the usual polar and axial angles with respect to the fixed {x,y, z} coordinate frame. In all simulations, 
we found very small changes in the magnitude of the individual spins up to the merger, hence these sky-map plots 
provide a very good representation of the spin dynamics. The left plot in Fig. [5] shows the individual spins, and 
the right plot shows the evolution of S. All the cases start with (pi = —180° and 1^92 = 0°. There are a couple of 
interesting aspects to notice in Fig. [8] First, there is no significant change in the ^1,2 direction, and hence no change 
in the 9 direction of S. Second, in all cases in the B-series, the precession is Ai^ « 120°. Since all the models start 
with the same 1^1,2, the spin orientation of the BHs arrive at the plunge (the point beyond which most of the kick is 
accumulated) with the same (p entrance angle. As we shall see in Sec. IVIl these two facts, particular to the B-series, 
have an important implication when fitting the gravitational recoils to the Kidder kick formula. 

Fig. [9] shows representative evolution tracks of the Si, 2 and S direction in the 9-ip plane for the S-series. The left 
and central plots in Fig. [Slshow the tracks of Si. 2 for some of the cases. The left plot includes the 0° < 6* = 6*2 < 180° 
models, with the central plot showing the 180° < 9 = 92 < 360° cases. The right plot in Fig. [S] depicts the evolution 
of S for the cases in the left and central plots. All the cases starts out ipi — —180°, (p2 = 0° and 9i = —90°. It is 
clear from Fig. [5] that the spin dynamics are significantly more complicated than in the B-series case. A substantial 
evolution in the 9 direction is evident in all cases, and there is also appreciable variation on the rate of (p precession 
from case to case. There is however a hint of a pattern. The closer the spin of BH2 aligns or anti-aligns with the 
z-axis, the larger is the evolution in the 9 direction. 
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FIG. 4: Same comparison as in Fig. |3] but for the model S-45. 




FIG. 5: Same comparison as in Fig. |3] but for the model S-90. 



VI. RECOIL ESTIMATES 

The gra vitational recoil from spinning BHs has been studied for a number of different initial spin configurations [isl . 
[TgI . ITtI . 113, IS] including very generic configurations [13] and for a systematic study of variations of the entrance 
angle in the orbital plane (i.e. xy-plane) between the spin vector and the a;-axis of anti-aligned BHs in Ref. [2Q|]. Our 
study explores the recoil of spin orientations out of the ccy-plane. Among other things, our aim is to test the assumption 
implied by the Kidder kick formula Eq. ^ that the recoil velocity can be split into components perpendicular and 
parallel to the orbital plane that depend on spin entrance angles at the plunge. 
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FIG. 6: Same comparison as in Fig. |3]but for the model S-135. 
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FIG. 7: Same comparison as in Fig. |3]but for the model S-225. 

We now specialize the Kidder kick formula to the B-Series. We denote by 9 the angle between a and the orbital 
angular momentum direction 1. In addition, the angle (p is the axial angle in the n-k plane relative to the n direction. 
In terms of these angles, the cosine directions in the Kidder kick formula Eq. ([5]) read: 



(la) 
(na) 



— cos 9 

= sin 9 cos (f 

= sin 9 sin . 
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FIG. 8: The evolution tracks of the Si, 2 and S directions in the d-ip plane for all the cases in the B-series. The left plot shows 
the individual spins and the right plot shows the evolution of S. All the cases start with 951 = —180° and ip2 — 0° . Notice that 
there is almost no change in the Oi^2 direction for the individual spins or for S. 




FIG. 9: Representative evolution tracks of the Si, 2 and S direction in the 9-ip plane for the S-series. The left plot shows the 
tracks of Si, 2 for some of the cases in which 0° < 9 — 62 < 180°, the central plot for 180° < 6 = 62 < 360° and the right plot 
the evolution of S for the cases in the left and central plots. All the cases start with i^i = —180°, 1^2 = 0° and 6\ — —90°. 



For generic cases, the angles 9 and cp are different from the polar angle 9 and axial angle (p introduced in Sec. [Til 
which were defined with respect to the fixed {x,y,z} coordinate system. This is because the {l,n, k} system, by 
design, is attached to the orbital motion of the binary; hence, it will follow also its precession. However, for all the 
cases we have considered, to a good approximation, the vector 1 stays aligned with the z-axis. Thus, 9 ~ 9. 

One of the goals of our work is to single out and explore the 9 projection dependence. That was the main motivation 
for constructing the B-series. In Sec. |Vl we saw that for each model in the B-series the angles ^1,2 remained fairly 
constant and the precession was such that the angles (^1,2 changed by the same amount in all models. As a consequence, 
it is possible to use in the Kidder kick formula and write the kick velocity in terms of the x, y and ^-components 
as: 

V = Co cos 9 
yy = Co Hy cos 9 

= CoK^s\^9 (16) 

where Co — Sg^/ (M^(l + g)^) and = Kk sin + Kn cos ip. and Hy are related to and Hk by a rotation in 
the xy-plane. Since in our study g = 1 and S/M^ 0.6, then Co 0.0375. Notice that V^,^^ = ^{9 = 0") = Co Hk, 
= VHO = 0°) = Co and V^^, ^V'{9 = 90°) = CoK. 
Fig. [TO] shows the x, y and z-components of the recoil velocity as a function of the initial value of 9 for all the cases 
in the B-series. We have also added the 9 = 0° case studied in Ref . [13] ■ The gravitational recoil was computed from 
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(degrees) 

FIG. 10: Recoil velocity components for the B-series as a function of the initial angle 6, the angle between S and the orbital 
angular momentum. Circles denote numerical values from the simulations, and the lines are obtained from l/^^'" ^ = V^^x ' cos f 
and ViJax sin(6'), where Vmajf^ are simply the recoil velocity components obtained for the B-0 and, similarly, Vj^ax for the B-90 
model. 

the Newman-Penrose quantity ^"4 at r/M — {30,40,50,60}. The plot shows r ~ 30 M. The results for the other 
detectors are of similar quahty except for r — 60 M where the resolution drops. In addition to the recoil data, we also 
shows the ^ cos 9 and V^l^^ sin 6 where Vi4ax ^ are simply the recoil velocity components obtained 

for the B-0 and, similarly, l^nax f^^' the B-90 model. We emphasize that no fitting to a s'm9 or cos 9 function was 
done in constructing Fig. (TUl Clearly the recoil velocity follows the sin 9 and cos 9 curves which is expected from the 
recoil formulas Eqs. (|16p and (jlip . This was possible because for the B-series there is a clear way of measuring the 
entrance angles. We found that V^^^ = 80 ± 12, V^^^ = 275 ± 41 and V^^^^^ = 960 ± 144 km s"^ , which yields the 
constants = (2.1 ± 0.3) • lO'', Hy = (7.3 ± 1) • lO'' and = (2.6 ± 0.4) • 10^ 

Because of the complicated dynamics in the S-series, we were not able to find a simple method for determining the 
entrance angles. As a consequence, it was not possible to do fittings to the Kidder kick formula. We are currently 
investigating [49j an approach that explicitly accounts for the precession dynamics that could potentially handling 
arbitrary configurations. 

VII. CONCLUSIONS 

The dynamics of BHs in interaction and merger, the gravitational radiation produced and the resulting kick in 
the final merge BH have direct implementations for understanding a wide range of astrophysical phenomena. These 
include the development of large scale structure, the structural evolution of galaxies, the detectability (for instance 
in the detector LISA) of gravitational radiation from the merger, and the statistics of double-nucleus galaxies. 

Our work concentrated on investigating the dynamics of spins in BBH systems and the gravitational recoil that 
the final BH experiences as a result of the merger. Regarding the spin dynamics, we have shown that spin precession 
follows fairly well the 2PN predictions up to the merger. Although we have only investigated two families of initial 
orientations {B-series and S-series), we believe that they represent fairly generic orientations, thus we speculate that 
the spin dynamics agreement with 2PN will be true for all orientation cases. It remains to be seen whether the 
agreement deteriorates when relaxing the condition of equal spin-magnitudes and/or masses. Spin-spin PN effects 
were not found to be significant for the cases we considered. 

An interesting aspect of the B-series, with BH spins initially anti-aligned with respect to each other, was that for 
each case the spins precessed about the orbital angular momentum axis, while keeping their polar (9) angle very 
closely constant. Also very interesting is that for all the models in the B-series, the vector S precessed almost the 
same amount about the orbital angular momentum axis. We were therefore able to read off the entrance angles and 
to demonstrate that the sin 9 and cos 9 dependences in the rate of linear momentum radiated as derived in Eq. (jlip 
get directly translated into the Kidder kick formula Eq. ([5]). 

For the S'-series a more complicated spin dynamics is found and the lack of symmetry between the BHs allows more 
complicated radiation and kick results in this case. We will continue addressing comparisons to 2PN and the validity 
of the Kidder kick formula for generic configurations in a separate paper [4^ . 
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